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We present a sharp-interface model of two-dimensional ramified growth during quasi-steady elec¬ 
trodeposition. Our model differs from previous modeling methods in that it includes the important 
effects of extended space-charge regions and nonlinear electrode reactions. The model is validated 
by comparing its behavior in the initial stage with the predictions of a linear stability analysis. 


I. INTRODUCTION 

Electrodeposition is a technologically important pro¬ 
cess with diverse applications and implications, e.g. for 
battery technology, electroplating, and production of 
metal powders and microstructures M. For well over 
a century it has, however, been known that the layer de¬ 
posited during electrodeposition is prone to morpholog¬ 
ical instabilities, leading to ramified growth of the elec¬ 
trode surface. Over the years, a large number of ex¬ 
perimental, theoretical, and numerical studies have been 
devoted to increasing the understanding of this ramified 
growth regime [12l-ll9l]. Big contributions to our under¬ 
standing of the growth process have come from diffusion- 
limited aggregation (DLA) models [2(1 j2l| and, more re¬ 
cently, phase-field models similar to those which have 
successfully been applied to solidification problems [22l — 
[26j . However, while both of these approaches capture 
parts of the essential behaviour of ramified growth, they 
have some fundamental shortcomings when applied to 
the electrodeposition problem. 

The first of these shortcomings has to do with the 
ion transport in the system. Typically, the electrolyte 
contains a cation of the electrode metal which can both 
deposit on the electrodes and be emitted from the elec¬ 
trodes. The anion, on the other hand, is blocked by the 
electrodes. The electrodes thus act as ion-selective ele¬ 
ments, and for this reason the system exhibits concen¬ 
tration polarization when a voltage is applied. In 1967, 
Smyrl and Newman showed [27 ] that in systems exhibit¬ 
ing concentration polarization, the linear ambipolar dif¬ 
fusion equation breaks down when the applied voltage 
exceeds a few thermal voltages. At higher voltages a 
non-equilibrium extended space-charge region develops 
next to the cathode, causing the transport properties of 
the system to change dramatically. It seems apparent 
that this change in transport properties must also lead 
to a change in electrode growth behavior. Indeed, this 
poi nt was argued by Chazalviel already in his 1990 paper 
112 . Now, the issue with DLA and phase-field models is 
that neither of these methods account for non-zero space- 
charge densities. It is therefore only reasonable to apply 
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these methods in the linear regime, where the applied 
voltage is smaller than a few thermal voltages. 

The other shortcoming of DLA and phase-field meth¬ 
ods is their treatment of the electrode-electrolyte inter¬ 
face. It is well known in electrochemistry that electrode¬ 
position occurs with a certain reaction rate, which is 
dependent on the electrode overpotential and typically 
modelled using a Butler-Volmer type expression [2j|, (29| . 
Nevertheless, electrode reactions are neither included in 
DLA methods nor in phase-field methods. 

There have been attempts to include finite space- 
charge densities in phase-field models, but the resulting 
models are only practical for ID systems because they re¬ 
quire an extremely dense meshing of the computational 
domain [22i, [23| . Attempts at including electrode reac¬ 
tions suffer from similar problems, as the proposed mod¬ 
els are sensitive to the width of the interface region and 
to the interpolation function used in the interface region 

mill. 

To circumvent the shortcomings of the established 
models we pursue a different solution strategy in this 
paper. Rather than defining the interface via a smoothly 
varying time-dependent parameter as in the phase-field 
models, we employ a sharp-interface model, in which 
the interface is moved for each discrete time step. Us¬ 
ing a sharp-interface model has the distinct advantage 
that electrode reactions are easily implemented as bound¬ 
ary conditions. Likewise, it is fairly straight-forward 
to account for non-zero space charge densities in a 
sharp-interface model, see for instance our previous work 

Refs, mm- 

Like most previous models, our sharp-interface model 
of electrodeposition models the electrode growth in two 
dimensions. There have been some experiments in which 
ramified growth is confined to a single plane and is ef¬ 
fectually two dimensional [TtI. [jjl [jit . However, for most 
systems ramified growth occurs in all three dimensions. 
There will obviously be some discrepancy between our 
2D results and the 3D reality, but we are hopeful that 
our 2D model does in fact capture much of the essential 
behavior. 

At this stage, our sharp-interface model is only ap¬ 
plicable once the initial transients in the concentration 
distribution have died out. In its current form the model 
is therefore mainly suitable for small systems, in which 
the diffusive time scale is reasonably small. We aim at 
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FIG. 1. Sketch of the initial geometry of the system. Two 
co-planar metal electrodes of width W are placed a distance 
of 2 L apart. The gap between them is filled by an electrolyte 
with cation concentration c+ and anion concentration c_. A 
voltage difference of Vo is applied between the electrodes. 


removing this limitation in future work. 


II. MODEL SYSTEM 

The model system consists of two initially flat paral¬ 
lel metal electrodes of width W placed a distance of 2 L 
apart. In the space between the electrodes is a binary 
symmetric electrolyte of concentration cq, in which the 
cation is identical to the electrode material. The elec¬ 
trodes can thus act as both sources and sinks for the 
cation, whereas the anion can neither enter nor leave the 
system. A voltage difference Vo (in units of the thermal 
voltage Vt = k B T/e) is applied between the two elec¬ 
trodes, driving cations towards the top electrode and an¬ 
ions toward the bottom electrode. A sketch of the system 
is shown in Fig. [T] 

By depositing onto the top electrode we ensure that 
the ion concentration increases from top to bottom, so 
we do not have to take the possibility of gravitational 
convection into account. To limit the complexity of the 
treatment, we also disregard any electroosmotic motion, 
which may arise in the system. We note, however, that 
the sharp-interface model would be well suited to investi¬ 
gate the effects of electroosmosis, since the space charge 
density is an integral part of the model. 


III. SOLUTION METHOD 



FIG. 2. Sketch of the electrode growth. The electrode surface 
at time ti is indicated with a full line. In the time step ti+i—ti 
an amount of material A L is deposited on the electrode. On 
basis of the deposited material the geometry at time £i+i is 
created (indicated with a dashed line). 


way of continuing from the old solution of the transport- 
reaction problem. One way of getting around this issue 
is to separate the time scales in the problem. More to the 
point, we assume that the growth of the electrode hap¬ 
pens so slowly, compared to the transport time scales, 
that the transport problem always is in quasi steady- 
state. By treating the transport-reaction problem as be¬ 
ing in steady state in each time step, a solution can be 
computed without reference to solutions at previous time 
steps. 

Obviously, the quasi steady-state assumption is flawed 
in the initial time after a voltage is applied to the sys¬ 
tem, as the application of a voltage gives rise to some 
transients in the transport problem. However, after the 
initial transients have died out the assumption is quite 
reasonable, except for the case of very concentrated elec¬ 
trolytes. To see that, we consider the thickness AL of 
the electrode growth in a time interval Af, 

AL = a 3 AtJ + , (1) 


where a 3 is the volume of a metal atom in the solid phase 
and J + is the current density of metal ions entering the 
electrode. The current density is on the order of the 
limiting current 2cqD + /L : so the time scale associated 
with an electrode growth of AL is 

AL LAL /ox 
1 ~ ~ 2D+c 0 a 3 ' ( } 


On the other hand, the transport time scale associ¬ 
ated with the distance AL is 


.AL 

Miff 


AL 2 
2 D + ' 


( 3 ) 


The basic idea in our solution method is to solve the 
transport-reaction problem for each time step, and then 
use the calculated currents to find the amount of mate¬ 
rial deposited at the electrode. Based on this deposition 
rate the geometry is updated, and the transport-reaction 
problem is solved for a new time step, as illustrated in 
Fig. m 

The major difficulty in employing this method is that 
when the geometry is updated the computational do¬ 
main is also remeshed, so there is no straight-forward 


The ratio of the transport time scale to the growth time 
scale is thus 


.AL 

Miff 

At 


AL . 

—coa , 


( 4 ) 


which is indeed very much smaller than unity. 

As mentioned above, our model does not apply to the 
initial time after the voltage is applied. To estimate how 
this impacts our results, we make a comparison of the im¬ 
portant time scales. The time it takes for the transients 
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to die out is given by the diffusion time, 

L 2 


fL — 

r diff — 


2 D 


+ 


( 5 ) 


The growth rate of the most unstable harmonic per¬ 
turbation to the electrode surface we denote r max (see 
Ref. [36j]), and from this we obtain an instability time 
scale, 


tu 


1 


It is apparent that if 


^diff ^inst, 


( 6 ) 


( 7 ) 


then nothing interesting happens to the electrode surface 
in the time it takes the transients to disappear. In this 
case our quasi-steady approach is therefore justified. 

Even if t^ iS ti ns t our approach may be justified. If 
the total deposition time is much larger than tj iff , then 
what happens in the time before the transients die out 
is largely unimportant for the growth patterns observed 
in the end. Thus, though the quasi-steady assumption 
seems restrictive, it actually allows us to treat a fairly 
broad range of systems. 


IV. GOVERNING EQUATIONS 


A. Bulk equations 


The ion-current densities in the system are given as 


■71 D±coc± , 

= ln(c±) + z± 4 f>, 


( 8 a) 

( 8 b) 


where D± are the diffusivities of either ion, Cq is the ini¬ 
tial ion concentration, c± are the concentrations of either 
ion normalized by co, /i ± are the electrochemical poten¬ 
tials normalized by the thermal energy k B T , and (f> is the 
electrostatic potential normalized by the thermal voltage 
Vt = k B T/e. In steady state the Nernst-Planck equa¬ 
tions take the form 


0 = -V- J±. 


(9) 


The electrostatic part of the problem is governed by the 
Poisson equation, 


2A = — p = —z+c + — Z-C-, 
where the Debye length A D is given as 


An — 


' k B Te^ 


2e 2 co 

At the electrodes the anion flux vanishes, 
n • J_ = 0, 


( 10 ) 


( 11 ) 


( 12 ) 


and the cation flux is given by a reaction expression 

n-J+ = -R. (13) 

Rather than explicitly modelling the quasi-equilibrium 
Debye layers at the electrodes, we follow Ref. [32| and 
implement a condition of vanishing cation gradient at 
the cathode, 


n ■ Vc_|_ = 0. 


(14) 


The last degree of freedom is removed by requiring global 
conservation of anions, 


I (c_ - l) dV = 0. 

Ja 


B. Reaction expression 


(15) 


We model the reaction rate using the standard Butler- 
Volrner expression 0 , 

R = ko _ e -7«-(l-a0 Ztf+V) ; ( 16 ) 

where ko is the rate constant of the reaction, V is the 
non-dimensionalized electrode potential, k is the surface 
curvature, a is the charge-transfer coefficient, and 7 is 
given in terms of the surface energy 7 , 


7 = 


a 3 7 

]~f- 


(17) 


Here, a 3 is the volume occupied by one atom in the solid 
phase. 7 K is thus a measure of the energy per atom 
relative to the thermal energy. 


V. NUMERICAL STABILITY 

Due to the surface energy term in the reaction expres¬ 
sion, the surface is prone to numerical instability. In 
an attempt to reach the energetically favorable surface 
shape, the solver will sequentially overshoot and under¬ 
shoot the correct solution. The fundamental issue we are 
facing is that the problem at hand is numerically stiff. As 
long as we are using an explicit time-integration method 
we are therefore likely to encounter numerical instabili¬ 
ties. 

The straight-forward way of updating the position r of 
the interface is to use the explicit Euler method, 


r(t + At) = r(t ) + na 3 AtR(t ), 


(18) 


where R(t) is the (position dependent) reaction rate at 
time t. To avoid numerical instabilities, we should in¬ 
stead use the implicit Euler method, 


r(t + At) = r(t) + na 3 AtR(t + At), 


(19) 
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where the reaction rate is evaluated at the endpoint in¬ 
stead of at the initial point. This is however easier said 
than done. R(t + At) depends on r(t + At) as well as 
on the concentration and potential distribution at t + At. 
Even worse, through the curvature R(t+At) also depends 
on the spatial derivatives of r(t + At). 

The way forward is to exploit that only part of the 
physics give rise to numerical instabilities. It is therefore 
sufficient to evaluate the problematic surface energy at 
t + At and evaluate the remaining terms at t. For our 
purposes we can therefore make the approximation 


R(t + At) « R(t : n(t + At)), (20) 

where k is the curvature. This does still make for a quite 
complicated nonlinear PDE, but we are getting closer to 
something tractable. The difference in curvature between 
t and t +At is small (otherwise we are taking too big time 
steps), so we can approximate 

R(t , k (t + At)) ss i?(t, k(£)) + R!(t, n(t)) Ak, (21) 

where R' denotes R differentiated with respect to k and 
An = n(t + At) — k( t). The curvature can be written as 


k = 


36 
ds ’ 


( 22 ) 


where 6 is the tangential angle of the interface and s is 
the arc length along the interface. We therefore have 

Ak = k(£ + At) — n(t) = (23) 

OS2 OS\ 


where we have adopted the shorthand notation 1 and 2 
for time t and t+At, respectively. The arc lengths si and 
S 2 will obviously differ for any nonzero displacement, but 
this is a small effect compared to the angle difference. As 
an approximation we therefore use S 2 ~ si and obtain 


A/- 


3(02-0i) 

ds 1 


(24) 


The tangential angle is a function of the surface 
parametrization, 


tan(0i) = 

OX 1 

For small displacements we can approximate 
dy 2 3{yi + Ay) 


tan(0 2 ) = 


3x2 9(* 1 + Ax) 


—!=f)' 

3 Ay 


3 Ay \ 


dy 1 
3x 1 


3x 1 


dyi _ 

3x 1 ^ 3x\ J 

dyi 3 Ax 
dxi 3x 1 


.. , 3Ay 3Ax 

= tan(0i) + 1 -tan(0i)- 


3x 1 


dxi 


(25) 


( 26 ) 


The difference in tangential angles can then be written 


02 — 0i = arctan 


tan(0i) 


3Ay 
3x 1 


tan(0i) 


3 Ax 
3x 1 


1 


1 + tan 2 (0i) 


3 Ay 
3x 1 


tan(0i) 


3 Ax 
3x 1 


-0i 

(27) 


Returning to the implicit Euler method Eq. m , we 
project it onto the normal vector to obtain 


A L = a 3 AtR(t + At) 

~ a 3 At [i?(t, k(£)) + R'(t, «(£)) Ak] , (28) 


where AL = n ■ [r(t + At) — r(t )]. The increments in 
the x and y directions are related to AL via 


Aa; = n x AL 1 Ay = n y AL. (29) 

Inserting these in Eq. (12711 and writing out the curvature 
difference Ak, we obtain a linear PDE for the displace¬ 
ment AL 


AL — a 3 AtR(t, n(t )) 
a 3 AtR’(t, k(£)) 

3 ( n y — n x tan(0i) 3 AL I 
3s\ \ l+tan 2 (0i) / 


(30) 


In the limit Ak = 0 this equation reduces to the original 
forward Euler method m- 


A. Correction for the curvature 


In the previous derivation, we did not take into account 
that the local curvature slightly changes the relation be¬ 
tween amount of deposited material and surface displace¬ 
ment AL. The deposited area in an angle segment d.0 can 
be calculated as 



71 


— T AL 

K 



d6 

Y 


AL 2 


AL 


(31) 


The line segment ds is related to the angle segment as 
ds = dd/ k. This means that 

o . . . gL4 k 

a 3 A tR(t + At) = — = - 
as 2 

= A L + ^A L 2 . (32) 

Using this expression in Eq. (l28l) yields the slightly non¬ 
linear PDE, with the term f kAL 2 , 

AL + (|AL 2 — a 3 AtR(t, k(£)) 
a 3 AtR' (t, k(£)) 

3 j n v — n x tan(0i) 9AL1 
9si \ l+tan 2 (0 i) 9*i J 

in place of Eq. 0 . 


AL 2 


AL 


( 33 ) 
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FIG. 3. (Color online) Three-dimensional extension of our 
two-dimensional model. The electrode interface can vary in 
the *y-plane according to the calculated ion-currents, but it 
has a fixed depth Ah in the 2 -direction. The interface is also 
divided into a number of bins of width As in the xy-plane. 
Each bin thus has the area AhAs. 


VI. NOISE 

An important part of the problem is the noise in the 
system, since the noise is what triggers the morphological 
instability and leads to formation of dendrites. Exactly 
how the noise should be defined is however a matter of 
some uncertainty. Most previous work uses a thermal 
white noise term with a small, but seemingly arbitrary 
amplitude. In this work we use a slightly different ap¬ 
proach, in which we assume that the noise is entirely 
attributed to shot noise. 

As it turns out, this approach requires us to be more 
specific about how our 2D model is related to the three- 
dimensional reality. In Fig. [3] a sketch of the tree¬ 
dimensional electrode is shown. The electrode interface is 
free to vary in the xy-plane, but has a fixed depth Ah in 
the ^-direction. Obviously, most real electrodeposits will 
have a more complicated behavior in the ^-direction, but 
for electrodeposits grown in a planar confined geometry 
this is actually a reasonable description. 

Solving the transport-reaction problem yields the cur¬ 
rent density at each point along the electrode surface, 
that is the average number of ions arriving per surface 
area per time. The mean number Q of ions arriving in 
an electrode section of size AhAs in a time interval At 
is thus 


Q = J+AhAsAt. (34) 

Since the ions are discrete entities, the actual number of 
arriving ions will, however, fluctuate randomly around 
the mean Q with some spread tx. We assume that within 
the time interval At, the arrival of each ion is statistically 
uncorrelated with the arrival of each other ion. It can 
then be shown that, as long as Q > 10, the number of 
arriving ions follow a normal distribution with mean Q 
and standard deviation 


ct = VQ. (35) 

This corresponds to an extra random current density 


■7rand 


VQ _ I J+ 
AhAsAt qrand ~ V A/iA.sAt ?rand ’ 


(36) 


where q ran d is a random number taken from a normal 
distribution with mean 0 and standard deviation 1. This 
in turn corresponds to a random electrode growth of 


A T _ 3.1 J+At _ 

^^rand — a \/ 


(37) 


Now, there is something slightly weird about this expres¬ 
sion for the random growth: it seems that the random 
growth becomes larger the smaller the bin size As is. 
However, as the bin size becomes smaller the weight of 
that bin in the overall behavior is also reduced. The net 
effect is that the bin size As does not matter for the 
random growth, see Appendix [A] for a more thorough 
treatment. 

The bin depth Ah, on the other hand, does matter for 
the random growth. Since our model is not concerned 
with what happens in the ^-direction, we simply have to 
choose a physically reasonable value of Ah, and accept 
that our choice will have some impact on the simulations. 
This is a price we pay for applying a 2D model to a 3D 
phenomenon. 


VII. NUMERICAL SOLUTION 

To solve the electrodeposition problem we use the 
commercially available finite element software COM- 
SOL Multiphysics ver. 4.3a together with MATLAB 
ver. 2013b. Following our previous work [H]', 32:, 37] , the 
governing equations and boundary conditions Eqs. (ISaf) . 
(HED, ©, m, (ECU m, m, ©, (USD, and © are 
rewritten in weak form and implemented in the mathe¬ 
matics module of COMSOL. For each time step the fol¬ 
lowing steps are carried out: First, a list of points defin¬ 
ing the current electrode surface is loaded into COM¬ 
SOL, and the surface is created using a cubic spline in¬ 
terpolation between the given points. The computational 
domain is meshed using a mesh size of As at the elec¬ 
trode surface, a mesh size of l in a small region next to 
the electrode, and a much coarser mesh in the remain¬ 
der of the domain. Next, the curvature of the surface 
is calculated at each point. The solution from the pre¬ 
vious time step is then interpolated onto the new grid, 
to provide a good initial guess for the transport-reaction 
problem. Then the transport-reaction problem is solved. 
Based on the solution to the transport-reaction problem 
the electrode growth A L is calculated by solving Eq. (l33l) 
on the electrode boundary. At each mesh point a small 
random contribution AL ran d = a 3 At J ran d is then added 
to A L. Finally, the new x and y positions are calculated 
by adding n x (AL + AL ran( j) and n y (AL + AL ran( j) to the 
old x and y positions. 

The new x and y positions are exported to MATLAB. 
In MATLAB any inconsistencies arising from the elec¬ 
trode growth are resolved. If, for instance, the electrode 
surface intersects on itself, the points closest to each other 
at the intersection position are merged and any interme¬ 
diate points are discarded. This corresponds to creating 
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FIG. 4. Example of the simplifying cutting procedure. The 
reduced interface (thick black line) divides the domain into 
an active region (white) and a passive region (light gray). 
The dark gray area shows the real cathode. The example is 
taken from a simulation with co = 1 mM and Vo = 10 after 
deposition for 31 hours and 28 minutes. 

a hollow region in the electrode which is no longer in 
contact with the remaining electrolyte. The points are 
then interpolated so that they are evenly spaced, and ex¬ 
ported to COMSOL so that the entire procedure can be 
repeated for a new time step. 

The simulations are run on a standard work station 
with two 2.67 GHz Intel Xeon processors and 48 GB 
RAM. The electrodeposits shown in Section IVIIII typi¬ 
cally take 2 days to run. 

A. Reduction of the computational domain 

At the cathode the mesh is much finer than in the 
remainder of the domain. The number of mesh points, 
and hence the computation time, therefore roughly scales 
with the length of the electrolyte-cathode interface. This 
has the unfortunate consequence that the computa¬ 
tion time for each time step increases drastically, when 
branching structures emerge at the cathode. To lower the 
computation time we exploit the fact that the vast ma¬ 
jority of the current enters near the tips of the dendritic 
structures. The parts of the cathode which are not near 
the tips can therefore be left fixed in time and thus re¬ 
moved from the simulation, without changing the results 
appreciably. This part of the domain is denoted the pas¬ 
sive region. In regions where the current density is less 
than 0.001 times the maximum value, we thus substitute 
the real, ramified electrode with a smooth line connect¬ 
ing the parts of the electrode with larger currents. The 
procedure is carried out in such a way that the real elec¬ 
trode surface can always be recovered from the reduced 
surface. For a few select examples we have verified that 
the results are unchanged by this simplifying procedure. 


TABLE I. Fixed parameter values used in the simulations. 


Parameter 

Symbol 

Value 

Cation diffusivity [38] 

D+ 

0.714 x 10“ 9 m 2 /s 

Anion diffusivity [381 

D _ 

1.065 x 10~ 9 nr 2 /s 

Ion valence 

Z 

2 

Surface energy 

7 

1.85 J/m 2 

Temperature 

T 

300 K 

Permittivity of water 

Gw 

6.90 x 10 _10 F/nr 

Charge-transfer coefficient 

a 

0.5 

Reaction constant^ 

k 0 

9.4 x 10 19 m“ 2 s _1 

Diameter of a copper atomb 

a 

0.228 nm 


a Calculated using the exchange current /q = 30 A/m 2 from 
Ref. j.hl and ko = Ia/(Ze). 

b The cubic root of the volume per atom in solid copper ioT . 


In Fig. 2] is shown an example electrode surface together 
with the reduced surface. It is seen that the length of 
the electrolyte-cathode interface is heavily reduced by ex¬ 
cluding parts of the electrode from the computation. 


B. Parameter values 

To limit the parameter space we choose fixed, physi¬ 
cally reasonable values for the parameters listed in Ta¬ 
ble |U The values are chosen to correspond to copp er elec¬ 
trodes in a copper sulfate solution, see Ref. j36| for de¬ 
tails. 

In Ref. [13 we calculate the critical wavelength A c , 
i.e. the smallest unstable perturbation wavelength, for a 
range of parameters. We expect the critical wavelength 
to be the smallest feature in the problem, so we choose 
the mesh size accordingly. We set the mesh size at the 
electrode to As = 0.1A C , since our investigations, see 
Section IVH Cl show that this is a suitable resolution. We 
also require that the mesh size does not exceed 0.1 times 
the local radius of curvature. In the bulk part of the sys¬ 
tem we use a relatively coarse triangular mesh with mesh 
size W/6. Close to the cathode, in a region l = 0.5 pm 
from the electrode surface, we use a triangular mesh with 
mesh size 1/4. See Fig. [5] for a meshing example. 

We choose a fixed value for the bin depth Ah = 0.2A C . 
In accordance with the analysis in Appendix [A] the time 
step At is chosen so that it is always smaller than 
0.5/r max . In addition, the time step is chosen so that 
at each point on the cathode, the growth during the time 
step is smaller than the local radius of curvature. 

We fix the length L to 100 pm. According to the time- 
scale analysis in Section mu and the instability growth 
rates found in Ref. 36j], the quasi-steady state approx¬ 
imation is valid for L = 100 pm. The width W of the 
system is set to W = 200A c , rounded to the nearest mi¬ 
crometer. This makes for a system that is broad enough 
to exhibit interesting growth patterns, while having a rea¬ 
sonable computation time. The growth is somewhat af¬ 
fected by the symmetry boundaries at y = 0 and y = W, 
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y L“ m l v [Mm] 

FIG. 5. (Color online) Example of domain meshing at varying 
magnification. The example is taken from a simulation with 
Co = 1 mM and Vo = 10 after deposition for 7 hours and 
50 minutes. The wiggly black line is the cathode surface. 
The light gray lines are the mesh boundaries and the dark 
(red) lines show the sections that are magnified. The mesh 
elements above the cathode surface are only used for storing 
the solution between time steps. 


especially at later times. 

These choices leave us with two free parameters, 
which are the bias voltage Vo and the electrolyte 
concentration Co- We solve the system for Co = 
{1 mM, 10 mM, 100 mM} and V 0 = {10,20,30}. 


C. Validation 


The random nature of the phenomena we are investi¬ 
gating poses obvious challenges when it comes to vali¬ 
dating the numerical simulations. The individual steps 
in the computation can be, and have been, thoroughly 
tested and validated, but testing whether the aggregate 
behavior after many time steps is correct is a much taller 
order. At some level, we simply have to trust that, if the 
individual steps are working correctly, then the aggregate 
behavior is also correct. To support this view, there is 
one test we can make of the aggregate behavior in the 
very earliest part of the simulation. 

In the early stages of the simulation the electrode sur¬ 
face is deformed so little, that the linear stability analysis 



FIG. 6. (Color online) Power spectra averaged over 50 runs 
for three different mesh sizes, As = {0.1A C , 0.2A C , 0.4A C }. In 
each run we used M = 100 time steps of At = 0.64 s and the 
parameter values Co = 10 mM, L = 100 pm, and Vo = 30. 
The full black line shows the analytical result and the dashed 
black lines show the analytical standard error on the mean. 
The result for As = 0.1 A c is shown in dark (red), the result 
for As = 0.2A C is shown in medium (red), and the result for 
As = 0.4A C is shown in bright (red). 


from [36| should still be valid. We thus have an analyti¬ 
cal expression for the wavelength dependent growth rate 
T, which we can compare with the growth rates found 
in the numerical simulations. In Appendix O we calcu¬ 
late an expression for the average power spectrum of the 
cathode interface after deposition for a time itot> given 
the type of noise described in Section IVTl 


(Pn) = 


J 4 


2AhWT ri 


O 2r n t tot 


- 1 ]. 


(38) 


where T n is the growth rate of the n’th wavelength \ n = 
W/n component in the noise spectrum. We also find the 
standard deviation SD(P n ) of the power spectrum 


SD(P„) « V2(P n ). (39) 


Because the standard deviation of P n is so large com¬ 
pared to the mean value, it is necessary to average over 
many runs before a meaningful comparison with Eq. (13811 
can be made. Averaging the power spectrum over 50 sim¬ 
ulations brings the standard error on the mean down to 
20 percent times the mean value, at which point a rea¬ 
sonable comparison can be made. In Fig. [ 6 ] the power 
spectrum averaged over 50 runs is shown for three dif¬ 
ferent mesh sizes, As = {0.1A C , 0.2A C , 0.4A C }. In each 
run we used M = 100 time steps of At = 0.64 s and 
the parameter values Co = 10 mM, L = 100 pm, and 
Vo = 30. The chosen step size corresponds to 0.01/T max . 
The analytical result (1381) is also shown together with 
the standard error on the mean. The power spectra are 
normalized with the power P° obtained for T = 0, 


p o _ „6 J+ttot 

AhlV' 


(40) 


It is seen that for As = 0.4A C some of the power in 
the small wavelength components is filtered out. As the 



























mesh size is decreased to As = 0.2A C and As = 0.1 A c the 
low wavelength components are represented increasingly 
well. 

In the above treatment, the time step was chosen 
very small compared to the instability time scale, At = 
0.01/r max . This was done to approach the limit of con¬ 
tinuous time, and thus enable the best possible compar¬ 
ison with the analytical theory. Such a short time step 
is, however, impractical for the much longer simulations 
in the remainder of the paper. In those simulations we 
use time steps as large as At = 0.5/r max . Due to the 
coarser time resolution employed in the remaining simu¬ 
lations, we expect their power spectrum to deviate some¬ 
what from the almost ideal behavior seen in Fig. [G] 


VIII. RESULTS 

We let the simulations run until the cathode has grown 
25 pm. The time to it takes to reach this point varies 
greatly with the parameters, mainly because the limiting 
current scales with Co- In Fig. [3 the cathode surfaces are 
shown along with heat plots showing the relative mag¬ 
nitude of the current density at the last time step. The 
white line shows the position of the reduced interface 
at the last time step, and the gray area shows the ac¬ 
tual position and shape of the cathode. The gray elec¬ 
trodeposits have different shades corresponding to 0.25fo, 
0.5fo, 0.75to, and t 0 . The heat plot shows the value of 
jn° rm , w } 1 j c } 1 j s the magnitude of the cation current den¬ 
sity normalized with its maximum value. In each panel 
jnorm t hus varies from 0 to 1. 

To investigate the reproducibility of the results we have 
repeated the simulation of the cq = 1 mM, Vo = 10 
system two times. All three electrodeposits are seen in 
Fig. [S] The electrodeposits are clearly different from one 
another, as expected for a random process, but they are 
also seen to share some general features. These shared 
features are most easily appreciated by comparing the 
electrodeposits in Fig. [8] to the electrodeposits in Fig. [3 
It is seen that the electrodeposits in Fig.[S]are much more 
similar to each other, than to any of the remaining elec¬ 
trodeposits in Fig. 0 Thus, the results are reproducible 
in the sense, that the random electrodeposits have some 
general features that are determined by the parameter 
values. 

When interpreting the plots in Fig. [3 we should be 
mindful that the aspect ratio is not the same in each 
panel. The reason for this is that the vertical axis has 
the same length, 30 pm, in each panel, while the length of 
the horizontal axis, W, varies between panels. In Fig. [9] 
we show adapted versions of the panels from Fig. [3 The 
subfigures in Fig. [9] are created by repeatedly mirroring 
the subfigures from Fig. [3 until their horizontal length is 
100 pm. Obviously, the resulting extended cathodes are 
somewhat artificial, since we have imposed some sym¬ 
metries, which would not be present in a simulation of 
a system with W = 100 pm. Nevertheless, we find the 


subfigures in Fig. [9] useful, since they give a rough im¬ 
pression of the appearance of wider systems and allow 
for easier comparison of length scales between panels. 


A. Rationalizing the cathode morphologies 

The cathode morphologies observed in Fig.[3and Fig. [9] 
are a function of several factors, some of which we at¬ 
tempt to outline below. First, we consider the time to 
it takes before part of the cathode reaches x = 175 pm. 
As seen from Eq. m, this time is mainly a function of 
the limiting current. This explains the approximately in¬ 
verse scaling with cq. The current density also increases 
with Vo, which is why the time to decreases slightly as 
Vo increases. Finally, the time to scales with the filling 
factor. This is the reason why to is much larger in the 
upper left panel of Fig. [3 than in either of the two other 
top row panels. 

It is apparent from the lack of ramified growth, that 
the cathode in the upper left panel in Fig. [3 is consider¬ 
ably more stable than the other systems in the leftmost 
column. To explain this variation in stability, we refer to 
Fig. 6 in Ref. [36|. There it is shown that the instability 
length scale is on the order of 50 pm for Co = 100 mM at 
Vo = 10, while it is considerably lower for Co = 10 mM 
and co = 1 mM. Fig. 6 in Ref. [36j also shows that for 
Vo > 18 the instability length scale decreases in size as 
the concentration increases. The same tendency is ob¬ 
served in Fig. [9] 

From the subfigures in Fig. [9] it appears that there 
is a connection between the thickness of the layer de¬ 
posited before the instabilities develop, and the charac¬ 
teristic length scale of the ramified electrodeposits. The 
analysis in Ref. [36j suggests that there is indeed such a 
connection and, moreover, that both lengths should scale 
with the most unstable wavelength for the given param¬ 
eters. To test this assertion, we plot the thickness <5i nst of 
the layer deposited before the instabilities develop, ver¬ 
sus the most unstable wavelength A max . We exclude the 
Co = 100 mM, Vo = 10 system, since instabilities have 
not yet developed in this system. The resulting plot is 
seen in Fig. EDA) together with a linear fit. Although 
there is a good amount of scatter around the linear fit, it 
is seen to capture the general trend reasonably well. 

We would like to make a similar plot with the charac¬ 
teristic length scale <5 c har of the ramified electrodeposits 
on the p-axis. To extract (5 c har, we follow the approach 
in Ref. [4C{ and calculate the so-called Minkowski dimen¬ 
sion of each electrodeposit. In doing this we only consider 
the part of the electrodeposit lying between 170 pm and 
190 pm, and as before we exclude the cq = 100 mM, 
Vo = 10 system. In this work we are actually not inter¬ 
ested in the Minkowski dimension itself, but rather in a 
partial result that follows from the analysis. In a range 
of length scales the electrodeposits appear roughly frac¬ 
tal, but below a certain length scale the electrodeposits 
are locally smooth. The length scale at which this tran- 
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FIG. 7. (Color online) Electrodeposits in the Vo-co plane obtained for L = 100 pm, Co = {1 mM, 10 mM, 100 mM) and 
Vo = {10,20,30}. The aspect ratio varies between the panels, since the width W of the simulated region is always set to 
200Amax- The gray area has different shades corresponding to times to (light), 0.75to (darker), 0.5to (darker yet), 0.25to 
(darkest). The white line indicates the reduced surface at time to- The contours in the liquid represent the relative magnitude 
of the cation current. 



170--- 

20 10 0 20 10 0 20 10 0 

y l/H y [aH y [/H 


FIG. 8. Three simulations of electrodeposits using the same 
parameter values L = 100 pm, Co = 1 mM, and Vo = 10. 
The electrodeposits are clearly different from one another, 
but they do share some general features. 


sition occurs can be extracted from the analysis, and we 
use this length as the characteristic length scale <5 c har of 


the electrodeposit, see Appendix iDl In Fig. fTOl b) we plot 
^char versus A max . Also here, we find a roughly linear 
behavior. 

Evidently, A max plays an important role for the mor¬ 
phology of the electrodeposits. However, <5j nst and <5 c har 
alone are not sufficient to characterize the electrode¬ 
posits. As seen in the top row of Fig. [9} the characteristic 
length scale <5 c har varies very little between Vo = 20 and 
Vo = 30. Yet, the morphology still changes appreciably. 
The reason for this change in morphology is probably 
that the gradient in electrochemical potential increases 
near the cathode as the bias voltage is increased. The 
larger the electrochemical gradient is, the more the sys¬ 
tem will favor deposition at the most protruding parts of 
the electrodeposits. For large voltages we therefore ex¬ 
pect long and narrow electrodeposits, whereas we expect 
dense branching electrodeposits for low voltages. 
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100 



10 Y 0 20 30 

FIG. 9. Extended electrodeposits in the Vo-co plane obtained by mirroring those from Fig. [7] in their symmetry axes until the 
width equals 100 pm. The dashed line indicates the first mirror plane, i.e. the part between the dashed line and y = 100 pm 
are obtained by repeating the part marked by a black line. 



100 200 300 400 500 600 

Amax [nm] 


FIG. 10. (Color online) (a): The instability length scale <$mst 
obtained from the simulations, plotted versus the most un¬ 
stable wavelength A max . Also, a linear fit highlighting the 
roughly linear dependence is shown, (b): The characteris¬ 
tic length scale 5 C h ar obtained from the simulations, plotted 
versus the most unstable wavelength A max . Also, a linear fit 
highlighting the roughly linear dependence is shown. 


IX. DISCUSSION 

Our model improves on existing models in three im¬ 
portant ways: it can treat systems at overlimiting cur¬ 
rent including the extended space-charge region, it allows 
for a proper reaction boundary condition, and it can be 
tested against results from sharp-interface stability anal¬ 
yses. Our model is, however, not without issues of its 
own. Perhaps the most apparent of these is the quasi- 
steady-state assumption. This assumption limits the ap¬ 
plicability of the model to short systems, in which the dif¬ 
fusion time is small compared to the deposition time, as 
discussed in Section Hill In principle the phase-field mod¬ 
els are superior to our model in this aspect, since they do 
not have this limitation. However, it is not of practical 
relevance, as all of the published phase-field simulations 
are for systems so short that the quasi-steady-state as¬ 
sumption is valid anyway f2~dl l2?ii |. 

It is well known, that the strong electric fields at the 
dendrite tips give rise to electroosmotic velocity fields 
in the system (4ll - l43| . To simplify the treatment and 
bring out the essential physics of electrodeposition, we 
have chosen not to include fluid dynamics and advection 
in our model. However, it is straightforward to include 
these effects, see for instance our previous work [32] . 

One of the main advantages the sharp-interface model 
has over the phase-field models, is that it allows for the 
implementation of proper reaction boundary conditions. 
The standard Butler-Volmer model used in this paper 
is a first step towards realistic reaction boundary condi¬ 
tions. As elaborated by Bazant in Ref. [2l|, there are 
other reaction models, such as Marcus kinetics, which 
might better describe the electrode reactions. Also, the 
standard Butler-Volmer model has the contentious as- 
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sumption that the overpotential is the total potential 
drop over both the electrode-electrolyte interface and 
the Debye layer. A more realistic approach might be to 
model the Debye layer explicitly or include the Frumkin 
correction to the Butler-Volmer model (idj j. Further¬ 
more, a proper reaction expression should take the crystal 
structure of the material into account. There are simple 
ways of implementing crystal anisotropy in the surface 
tension term, see for instance Refs. [26j, HU, but again, 
to keep the model simple we have chosen not to include 
anisotropy at the present stage. Any of the above men¬ 
tioned reaction models can be easily implemented in the 
framework of the sharp-interface model, and as such the 
specific Butler-Volmer model used in this work does not 
constitute a fundamental limitation. 

More broadly, our sharp-interface model includes, or 
allows for the easy inclusion of, most effects that are im¬ 
portant for electrodeposition in 2D. A natural next step 
is therefore to see how our results compare to experi¬ 
mental electrodeposits. Unfortunately, most such exper¬ 
imental data are viewed at the millimeter or centimeter 
scale, whereas our simulation results are at the microm¬ 
eter scale. In one paper, Ref. jij], the electrodeposits are 
probed at the micrometer scale, but the results do not 
make for the best comparison, since the morphology of 
their electrodeposits was a result of adding a surface ac¬ 
tive molecule. We hope that as more experimental results 
become available, it will be possible to perform rigorous 
tests of our model. 


X. CONCLUSION 

We have developed a sharp-interface model of elec¬ 
trodeposition, which improves on existing models in a 
number of ways. Unlike earlier models, our model is able 
to handle sharp-interface boundary conditions, like the 
Butler-Volmer boundary condition, and it readily deals 
with regions with non-zero space-charge densities. A fur¬ 
ther advantage is that our model handles the physical 
problem in much the same way as done in various linear 
stability analyses. We can thus obtain a partial valida¬ 
tion of our model by comparing its predictions with those 
of a linear stability analysis. As of now, the main weak¬ 
ness of our model is that it assumes quasi-steady state 
in the transport equations. For the systems studied in 
this paper this is a reasonable assumption, since the dif¬ 
fusion time is small compared to the instability time. In 
future work we want to extend the model to the transient 
regime, so that larger systems can be treated as well. 

The main aim of this paper has been to establish 
the sharp-interface method, but we have also included 
a study of the simulated electrodeposits. An interesting 
observation is, that the characteristic length scale of the 
electrodeposits seems to vary linearly with the size of the 
most unstable wavelength. This exemplifies a promis¬ 
ing application of our sharp-interface model, namely as 
a tool to develop a more quantitative understanding of 


electrodeposits and their morphology. 


ACKNOWLEDGMENTS 

We thank Edwin Khoo and Prof. Martin Z. Bazant for 
valuable discussions of electrode reactions and the growth 
mechanisms. 


Appendix A: Initial growth 


In the initial part of the simulation the electrode is so 
flat that the linear stability analysis from Ref. (36| gives 
a good description of the growth. We parameterize the 
cathode position as 


x = X(t) + f(y,t), (Al) 

where f(y , t) is the y-dependent deviation from the mean 
electrode position X(t). According to the linear stability 
analysis each mode grows exponentially in time with the 
growth factor F. After a time t an initial perturbation, 

N 

f{y, 0) = 'Y, a n e iknV , (A2) 

n= 1 

has therefore evolved to 

N 

f(y, t) = Y a n e rnt e lknV . (A3) 

n= 1 


We note that some of the growth rates T n can be nega¬ 
tive. In our simulation we add new perturbations with 
small time intervals, which we, for the purpose of this 
analysis, assume to be evenly spaced. After M time in¬ 
tervals At the surface is therefore described by 

M N 

f(y,MAt) = —m)At gik n y ' (A4) 

m— 0 n= 1 


We are interested in the average power of each mode 


(Pn) = 


M 

771=0 


0 F n (M-m)At 


(A5) 


The coefficients are random and uncorrelated with zero 
mean. On average the cross-terms in the sum therefore 
cancel and we can simplify, 


M 


(p~) = E 




\ 771=0 


M 


= <K| 2 ) Y, e 2r ^ M ~ m ^ At 

771=0 

„2r„(M+l)At _ -1 

= (Kl 2 )- 


D 2r„ At 


- 1 


(A6) 
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The variance of the power is given as 

Var(P„) = <P 2 } - (P n > 2 . 
The first of these terms is 

M 

(Pi) = ( ( E 


(A7) 


= e 


m —0 


,4r n MA£ 


M 


^ ^ ^nmQ 
m =0 


(A8) 


where g = e r " At . Writing out the absolute value 


M M 


<^) = 


E E anma* nrn , q m+m ' 


\m '=0 m—0 


(A9) 

where superscript * denotes complex conjugation. Be¬ 
cause the coefficients are uncorrelated with mean 0, only 
the terms including \a nm \ 2 \a nm '\ 2 survive in the average 
of the square, 


(P 2 ) = e 4r " MAt / I E E Wi v (m+m,) ) 

\ rn'—O m—0 / 

M M 

= 3e 4r„MA t £ E <l a nm| 2 |anm'| 2 ) 9 2(m+Tn,) - 

m'=0 m—0 

(A10) 

Here, the factor of six comes from the binomial coefficient 
and the factor of a half takes into account that the double 
sum counts each combination twice. Now, there are two 
possibilities; either m / or m = m!. In the first case 
\a nm \ 2 and \a nm ’\ 2 are uncorrelated, meaning that 

(|<w| 2 |<W'| 2 ) = (|a n | 2 ) 2 ■ 

Whereas if m = m', then 


M M 


(\a n m\ 2 \a n m'\ 2 ) = (W 4 ) . 

This means that 

M M 

(P») = 3e«--»“'<|a n | 2 > 2 V Y.1 

m'^m m =0 
M 

+ 3 e^“*(\a n \*)Y l Q 


(All) 

(A12) 


t 2{rn-\-m) 


Am 


m—0 
M M 


= 3e 


4r«MAt 


(W ! )’EE« 


2(m+m / ) 


m'—O m =0 


T 3e 
= 3 (P n ) 2 


AF n MAt 


M 

(<K| 4 )-<K| 2 ) 2 ) E? 4m 


m=0 




2 \ g 4 (M+l) _ J 


q 4 - 1 


(A13) 


The variance of the power is thus given as 

2 . e 4r„(M+i)At _ 1 
) e 4r„At _ i 
(A14) 

If T n At < 1 we can expand the denominators of (P„) 2 
and the last term. We find that they scale as 4(T„Af) 2 
and 4T n Af, respectively. In the limit T n At -C 1 the 
first term thus dominates over the second, so to a good 
approximation we have 

Var(P„) « 2(P„) 2 , (A15) 

SD(P„) » y/2(P n ). (A16) 

In the simulations the surface perturbations have the 
form 


Var(P„) = 2(P„) 2 + (<|«„| 4 ) - (|a n | 2 ) 


where, 


N 

f(y, o) = E bnh (v ~ nA y)’ 

n= 1 



1, 0 < y < As, 
0 , else. 


(A17) 


(A18) 


We take the absolute square of f(y, 0) given as both 
Eq. (|A2ll and Eq. (1A17I) . and integrate over the domain 
to obtain 


pW a „w 

/ 1/(2/, 0)| 2 Ay = E Kf / \Hy ~ nAy)\ 2 dy 

J o n—l 


N 


r W 


n =1 

N r w 


= As E \bn\ 2 , 

n= 

N 

1 /(2/, 0) | 2 dy = E K 

n—l 

N 

= wj2 


(A19) 


2 / |gifc n y |2 


dy 


(A20) 


The mean square of b n is thus related to the mean square 
of a n as 


As 1 

(Kl 2 > = ^{\K\ 2 ) = ^( 1 b n \ 2 ). 


From Eq. (l37l) we have that 


/17 12 \ 6 J+^t 


Inserting in Eq. 


we find 


1 s J + At e 2r "( M+1 ) A ‘ - 1 

\Ml) — - - a 


N AhAs 


„ 2 r„At _ i 


6 J + At e 2r n( 4 tot+ A t) _ i 
° AhW i^r^A t~\ 


(A21) 


(A22) 


(A23) 

(A24) 
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V [aH 

FIG. 11. The box-counting method illustrated on the elec¬ 
trodeposit obtained for co = 10 mM and Vo = 20. The boxes 
that cover part of the deposit perimeter are shown in gray 
and the remaining boxes are shown in white. In this example 
the grid size is e = 0.85 pm and the number of boxes it takes 
to cover the perimeter is N(e) = 234. 


which is seen to be independent of the bin size As. We 
also introduced the total time f t ot = MAt. In a consis¬ 
tent scheme the power spectrum should of course only 
depend on the total time, and not on the size At of the 
time steps. For small values of r„At we can expand the 
denominator and neglect the At in the nominator, 


(Pn) 


J 4 




2AhWT ri 


-!]• 


(A25) 


So, as long as 2r„At -C 1 the power spectrum does not 
depend on the size of the time step. 

For larger values of 2r„At the power spectrum does 
depend on the size of the time step. However, as long 
as 2r„At < 1, we do not expect the overall morphology 
of the electrode to have a significant dependence on the 
time step. 



FIG. 12. The number N(e) of boxes it takes to cover the 
electrodeposit plotted vs the box side length e. A linear fit is 
shown in each of the two approximately linear regions, and the 
Minkowski dimension in each region is indicated. The crossing 
point between the linear fits is marked by an arrow, and the 
characteristic dimensions 5 C h ar = 0.50 pm is calculated based 
on this crossing point. 


For a proper fractal geometry, the Minkowski dimen¬ 
sion is defined as 


r ln[A(e)] 
hm 

—>-0 ln(e) 


(Bl) 


The electrodeposits we are investigating are not fractal 
at all length scales, but in a range of length scales, we 
can calculate an approximate Minkowski dimension as 
the negative slope in a In [7\T(e)] vs ln(e) plot. In Fig. [TS1 
such a plot is seen, together with linear fits in each of 
the two approximately linear regions. The Minkowski 
dimension at small e is nearly unity, indicating that the 
deposit perimeter is locally smooth at this length scale. 
For larger values of e the Minkowski dimension deviates 
from unity, because the deposit is approximately fractal 
in this size range. At the transition point between these 
two regions is the smallest length scale, which is related to 
the morphology of the electrodeposit. This length scale 
we denote the characteristic length <5 c har- Technically, we 
define <5 c har as the point where the linear fits from each 
region cross each other, as indicated in Fig. fH?l 


Appendix B: Characteristic length scale 

To find the characteristic length scale <5 c har of the ram¬ 
ified electrodeposits we follow Ref. [!o| and use the box¬ 
counting method to calculate the Minkowski dimension 
of the deposits. We place a square grid with side length 
e over each deposit, and count the number N(e) of boxes 
it takes to completely cover the perimeter of the part of 
the deposit lying between x = 170 pm and x = 190 pm. 
An example is shown in Fig. 1111 
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